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ABSTRACT 

The analytical theory of diffusive cosmic ray acceleration at parallel stationary shock waves with magneto¬ 
static turbulence is generalized to arbitrary shock speeds V s = Pic, including in particular relativistic speeds. 
This is achieved by applying the diffusion approximation to the relevant Fokker-Planck particle transport equa¬ 
tion formulated in the mixed comoving coordinate system. In this coordinate system the particle’s momentum 
coordinates p and p = p /p are taken in the rest frame of the streaming plasma, whereas the time and space co¬ 
ordinates are taken in the observer’s system. For magnetostatic slab turbulence the diffusion-convection trans¬ 
port equation for the isotropic (in the rest frame of the streaming plasma) part of the particle’s phase space den¬ 
sity is derived. For a step-wise shock velocity profile the steady-state diffusion-convection transport equation 
is solved. For a symmetric pitch-angle scattering Fokker-Planck coefficient D^—p) = D m (pj) the steady-state 
solution is independent of the microphysical scattering details. For nonrelativistic mono-momentum particle 
injection at the shock the differential number density of accelerated particles is a Lorentzian-type distribution 
function which at large momenta approaches a power law distribution function N(p > p,j °c p ’ with the spec¬ 
tral index ^((3 1 ) = 1 + [3/(Fi r 2 — pj — 1)(1 +3(3|)]. For nonrelativistic (pi <C 1) shock speeds this spectral 

index agrees with the known result q(p| <C 1) — (r + 2)/(r— 1), whereas for ultrarelativistic (Ti 2> 1) shock 
speeds the spectral index value is close to unity. 

Subject headings: acceleration of particles - cosmic rays - relativity - shock waves 


1. INTRODUCTION 

One of the most important problems of modern astro¬ 
physics is to explain how cosmic ray particles are acceler¬ 
ated to relativistic energies in powerful sources of nonther- 
mal radiation. Diffusive first-order Fermi acceleration at non¬ 
relativistic shock fronts has been regarded as a prime candi¬ 
date for particle acceleration in astrophysics (for reviews see 
Drury 1983; Blandford and Eichler 1987). Modern TeV air- 
Cherenkov telescopes have indeed resolved the shock regions 
in supernova remnants and identied the shocks as strong emis¬ 
sion regions of TeV photons generated by the accelerated par¬ 
ticles (Hinton and Hofmann 2009). 

It is likely that diffusive shock acceleration also operates ef¬ 
ficiently at magnetized shock waves with relativistic speeds. 
Such relativistic shocks form during the interaction of rela¬ 
tivistic supersonic and super-Alfvenic outflows with the am¬ 
bient ionized interstellar or intergalactic medium (Gerbig 
and Schlickeiser 2009) producing anisotropic counterstream 
plasma distribution functions due to shock-reflected charged 
particles in the upstream medium (Spitkovsky 2008, Sironi 
and Spitkovsky 2009). Relativistic outflows are a direct con¬ 
sequence of violent explosive events such as in gamma-ray 
burst sources (Piran 1999) both in the collapsar (Woosley 
1993, Paczynski 1998) and supranova (Vietri and Stella 1998) 
models, but also occur as highly collimated pulsar winds and 
jets of active galactic nuclei with initial bulk Lorentz factors 
r 0 = (1 - (Vb/c) 2 )" 1/2 ^ 400. 

The transport and acceleration of energetic particles in the 
partially turbulent cosmic magnetic fields associated with 
shocks is described using the Fokker-Planck equation for 
the particle distribution function (for a recent derivation, see 
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Schlickeiser 2011). The diffusion approximation for the par¬ 
ticle density in the rest frame of the fluid is a well-known sim- 
plied form of the Fokker-Planck equation, which results when 
turbulent pitch-angle scattering is strong enough to ensure 
that the scale of the particle density variation is signicantly 
greater than the particle mean free path (Jokipii 1966, Has- 
selmann and Wibberenz 1968, Earl 1974). Numerical studies 
confirmed the accuracy of the diffusion approximation in a 
uniform mean guide magnetic field (Kota et al. 1982). 

While for nonrelativistic shock waves the analytic theory 
of diffusive shock acceleration is well developed (Axford 
et al. 1977, Krymsky 1987, Blandford and Ostriker 1978, 
Bell 1978, Drury 1983), for relativistic shock speeds such 
an analytical theory does not exist sofar even for parallel 
shock waves, although the underlying Fokker-Planck trans¬ 
port equation (see Eq. below) for the particle dynam¬ 
ics has already been derived (Webb 1985; Kirk, Schneider 
and Schlickeiser 1988) some years ago. The existing litera¬ 
ture concentrated on semi-numerical eigenfunction solutions 
of this Fokker-Planck transport equation pioneered by Kirk 
and Schneider (1987), relativistic Monte Carlo simulations 
(Ellison et al. 1990, Ostrowski 1991, Bednarz and Ostrowski 
1998, Summerlin and Baring 2012), and relativistic particle 
in-cell simulations (Spitkovsky 2008, Sironi and Sptkovsky 
2009). It is the purpose of this manuscript to develop an an¬ 
alytical study of cosmic ray acceleration in parallel relativis¬ 
tic magnetized shock waves employing the diffusion approxi¬ 
mation in the upstream and downstream regions of the shock 
wave. The development runs much in parallel with the exist¬ 
ing work on nonrelativistic shocks. 

2. BASIC EQUATIONS 

Magnetized space plasmas such as the interstellar medium 
harbour low-frequency linear (8 B <C Bq) transverse MHD 
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waves (such as shear Alfven and magnetosonic plasma waves) 
with dispersion relations co| = V^k 2 and ©\ = Vyk 2 , respec¬ 
tively, in the rest frame of the moving plasma. Faraday s in¬ 
duction law then indicates for MHD waves that the strength of 
turbulent electric fields 8 E = {V.\ /c)8B -C 8 B is much smaller 
than the strength of turbulent magnetic fields. The ordering 
Bq ^ 8 fi 8 E corresponds to the derivation of cosmic ray 
transport equations for < /> (X,p,p,ty,t) —► fo(X,p,p,t) —>■ 
E(A,p,f) from the collisionfree Boltzmann equation fir the 
full phase space distribution </> (X,p,p,§,t) totheFokker- 
Planck equation for its gyrotropic part fo(X,p,p,t), and to the 
diffusion-convection transport equation for its isotropic part 
F(X,p,t), respectively (for a recent review see Schlickeiser 
2011). Accordingly, the cosmic ray anisotropy, defined as the 
deviation 


g(X,v,p,t)=fo(X,p,p,t) — F(X,p,t), (1) 


c 2 dF(z) dU ^3 d(UT) 

U (z) dz dz dz 


(4) 


In the Fokker-Planck Eq. (Rli the phase space coordinates 
have to be taken in the mixed comoving coordinate system 
(time and space coordinates in the laboratory (=observer) sys¬ 
tem and particle’s momentum coordinates p and p = p /p 
in the rest frame of the streaming plasma). Moreover, in 
Eq. 0 we use the Einstein sum convention for indices, 
and y ()) . y„ £ \p,X,Y] represent the three phase space vari¬ 
ables with non-vanishing stochastic fields 8 E and 8 B. Conse¬ 
quently, the term on the right-hand side generally represents 
9 different Fokker-Planck coefficients: but, depending on the 
turbulent fields considered, not all of them are non-zero and 
some are much larger than others. S(X,p,t) accounts for ad¬ 
ditional sources and sinks of particles. 

The focusing length (Roelof 1969) 


then is small (|g| <C F) with respect to F. The diffusion ap¬ 
proximation applied to the Fokker-Planck transport equation 
for /o(X,/?,/j,f) allows us to relate the cosmic ray anisotropy 
g to the solutions of the diffusion-convection tranport equa¬ 
tion for F. 

Because of the gyrorotation of the cosmic ray particles in 
the uniform magnetic field, one is not so much interested in 
their actual position as in the coordinates of the cosmic ray 
guiding center 


X = (X,T,Z) =r+ApxBo=r+4" ( ~Px ] , (2) 
q u B-Q q«B 0 \ o J 

where we orient the large-scale guide magnetic field, which 
is uniform on the scales of the cosmic ray particles gyradii 
Rl = v/|n|, Z?o = B()C Z = ( 0 , 0 , Bo) along the z-axis. v and 
FL a = q«Bo /ym a c denote the speed and the relativistic gyrofre- 
quency of a cosmic ray particle with mass m, charge g a and 
energy ym a c 2 . 

The Larmor-phase averaged Fokker-Planck transport equa¬ 
tion in a medium, propagating with the stationary bulk speed 
U = U(z)e z with F = [1 — (U 2 /c 2 )]~ 1 ^ 2 aligned along the 
magnetic field direction is given by (Webb 1985; Kirk, 
Schneider and Schlickeiser 1988) 
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irrespective of how the Fokker-Planck coefficients are calcu¬ 
lated, either by quasilinear (Schlickeiser 2002) or nonlinear 
(Shalchi 2009) cosmic ray transport theories. a(z) denotes the 
rate of adiabatic deceleration/acceleration in relativistic flows 


L= _^ _ B o _ -_ (5) 

Bq • grad (|B 0 |) ( dB 0 (z)/dz ) 

represents the spatial gradient of the guide magnetic field Bo¬ 
lt is the only term, resulting from the mirror force in the large 
scale inhomogeneous guide magnetic field, that we consider 
(neglecting possible drift effects). L > 0 for a diverging guide 
magnetic field; L < 0 for a converging guide magnetic field. 

= -P~ 2 ^ [P 2 Floss/] + y (6) 

represents continuous (/)| 0S s) and catastrophic (T c ) momentum 
losses of cosmic ray particles. 

For spatially constant flows the rate of adiabatic decelera¬ 
tion/acceleration |4} vanishes, and the remaining flow veloc¬ 
ity ((/) dependent terms in Eq. (Rl simply result from the 
Lorentz transformation of special relativity of the laboratory- 
frame position-time coordinates (z.- t ) to the mixed-comoving- 
frame position-time coordinates (z ,f ). However, for spa¬ 
tially varying flow speeds U (z) special relativity no longer ap¬ 
plies and has to be replaced by the transformation laws from 
general realtivity. As noted by Riffert (1986) as well as Kirk, 
Schneider and Schlickeiser (1988) these introduce connection 
coefficients or Christoffel symbols of the first kind. In a flat 
Euclidean space-time metric the terms proportional to a(z) in 
Eq. ([3]) are exactly these connection coefficients. 

3. MAGNETOSTATIC SLAB TURBULENCE 

As important special case we consider magnetostatic (van¬ 
ishing turbulent electric fields), isospectral, slab (k = 0) tur¬ 
bulence. Then the only nonvanishing Fokker-Planck coeffi¬ 
cient is D^. In this case the Fokker-Planck transport equation 
([3]> simplifies to 
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Due to the rapid pitch angle scattering the gyrotropic par¬ 
ticle distribution function fo(X,p,p,t) adjusts very quickly 
to a distribution function which is close to the isotropic 
distribution in the rest frame of the moving background 
plasma. Defining the isotropic part of the phase space den¬ 
sity F(X,z,p,t) as the //-averaged phase space density 


3.1. Diffusion approximation for cosmic shock waves 

With |g| <F we approximate the anisotropy equation © 
to leading order by 


Yvp 



U dF\ 
c 2 dt ) 



Up 
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dF 
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F(X,p,t) = dpf 0 {X,p,p,t), (8) 

we follow the analysis of Jokipii (1966) and Hasselmann and 
Wibberenz (1968) to split the total density /o into the isotropic 
part F and an anisotropic part g, 

f 0 (X,p,p,t) =F(X,p,t)+ g(X,p,p,t) (9) 

where because of Eq. 0 

f dpg(jt,p,p,t) = 0 (10) 

Inserting the ansatz d9ll yields for the Fokker-Planck equation 
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( 11 ) 


Averaging this equation over p, using that D llu (p = ±1) = 0, 
leads to the diffusion-convection transport equation 


d_ 
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dg 
D — 
m dp 


(14) 


In previous nonrelativistic (r ~ 1) studies (Jokipii 1966, Has¬ 
selmann and Wibberenz 1968, Schlickeiser et al. 2007, 
Schlickeiser and Shalchi 2008) only the first term on the left- 
hand side of Eq. © has been considered, so that 
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For gradual flows with non-zero values of a for an extended 
region of space, we have to consider also the term resulting 
from the momentum gradient {dF/dp). This case will be con¬ 
sidered elsewhere. 

Here we consider flows with non-zero values of a in a very 
limited region of space such as cosmic shock waves. We 
chose as laboratory frame the rest frame of the shock wave 
with the step-wise velocity profile 


U(z) = 


—U i = const, for 0 < z < °° (upstream) 

—U 2 = const, for — °° < z < 0 (downstream) 

(16) 

with Ui <U\. In this case the rate of adiabatic acceleration 
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( 12 ) 

involving the first and second moment of the anisotropy. 

Subtracting Eq. <n"2]) from the Fokker-Planck equation (fTT 1 
provides 


is non-zero only at the position of the shock. We therefore 
neglect in the anisotropy equation ( fl4| the term resulting from 
the momentum gradient (dF /dp), and additionally neglect the 
time derivate of F as compared to the spatial gradient of F, i.e. 


U dF dF 
c 2 dt ^ dz' 


(18) 


so that Eq. © reduces to the streaming anisotropy contribu¬ 
tion only 
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We note that Eqs. ([T2| and ([T3 |i are still exact. 



dpp 2 g 


(13) 


which differs from the nonrelativistic anisotropy equation 
(TTST) by the additional factor T. 

Integrating Eq. © provides 
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The integration constant co (with respect to p) is determined 
from the property that the left-hand side of this equation van¬ 
ishes for p = ±1, yielding cq = —(Yv/2)(dF/dz), providing 


dg Ev(l —p 2 ) dF 
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A further integration over p together with Eq. © provides 
for the cosmic ray anisotropy 
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which is determined by the spatial gradient of the isotropic 
distribution function F(X,Y,z,p,t ) with respect to z. 

3.2. Anisotropy moments 

With the anisotropy ( [22] ) we readily determine the moments 
needed in the diffusion-convection transport equation ( [T 2 ] i as 

TvK q dF 


with the cosmic ray particle Lorentz factor y = (1 p 2 ) 12 

and the dimensionless particle velocity p = v/c, the diffusion 
operator (28 i reads 
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With the moments pTj ) - ( | 22 ] > we derive after straightforward 
but tedious algebra 


-l 6 dz 

in terms of the two {n = 0 , 1 ) integrals 
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For the often considered case of symmetric pitch-angle 
Fokker-Planck coefficients D^(—p) = D^ffp) the integral 
= 0 vanishes. 

3.3. Diffusion-convection transport equation 

For further reduction of the diffusion-convection transport 
equation ( fT 2 | ) it is convenient to introduce as new momentum 
variable 

s = lnp (26) 

and to define the convection and diffusion operators 
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Likewise, the convection operator (27 1 becomes with the first 
moment 


rp = ric a dF a dF 
L dz 3 ds 
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Inserting the operators ( [32] ) and ( [34]) p rovides for the 
diffusion-convection transport equation (|29ji 
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Applying the identity 
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respectively. With the assumption ( [T8| we then can write the 
diffusion-convection transport equation ( 121 in the compact 
form 


for any momentum dependent quantity M(s) to the last term 
on the left-hand side of Eq. we obtain 
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Using for any momentum dependent quantity M(s) the iden¬ 
tity 


Consequently, we find as alternative form for the diffusion- 
convection transport equation ([35]) 
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With the diffusion coefficients ( |33[ i we obtain for the differ¬ 
ence 
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where we used Eq. 0. The diffusion-convection transport 
equation (381 then becomes 
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We note the identity 
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With these two identities the diffusion-convection transport 
equation (401 finally reads 
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+ 3LF = S(X,p,t ) (45) 


Eq. ( |45| > is the first important new result of this study: 
the diffusion-convection transport equation of cosmic rays in 
aligned parallel flows of arbitrary speed containing magneto¬ 
static slab turbulence with the cosmic ray phase space coor¬ 
dinates taken in the mixed comoving coordinate system. It 
is particularly appropriate to investigate cosmic ray particle 
acceleration in parallel relativistic flows. 

In the limit of nonrelativistic flows U (z) c so that T ~ 1. 
the transport equation (451 becomes 


dF K- z dF 


_LF_ 

p 2 dp 


dt 


P 2 K 


pz 


L dz 

dF 

dz 


d_ 

dz 


a p 2 F 


dF 

UF k zz — 


+ %F = S(X,p,t), (46) 


which differs from the transport theory used in earlier non¬ 
relativistic diffusive shock acceleration theory (Axford et al. 

1977, Krymsky 1987, Blandford and Ostriker 1978, Bell 

1978, Drury 1983) by the additional third last term on the 
left-hand side involving K pz , which results from our correct 
handling of the connection coefficients in Eq. 0. As we 
will demonstrate below, this additional term provides a ma¬ 
jor modification of the resulting differential momentum spec¬ 
trum of accelerated particles in the nonrelativistic flow limit 
at nonrelativistic particles momenta: instead of a power law 
distribution of accelerated particles at the shock a Lorentzian 
distribution function results, which at large momenta then ap¬ 
proaches the power law distribution inferred in earlier accel¬ 
eration theories for nonrelativistic shock speeds. 

3.4. Slab Alven waves 

There are four different magnetostatic slab Alven waves: 
forward and backward propagating waves which each can 
be left- or right-handed circularly polarized, respectively. In 
terms of the cross helicity H c and the magnetic helicities o± 
of forward and backward moving Alfven waves, and power- 
law type wave intensities I « fcjp with sS (1,2) with the same 
spectral index s of all four waves (isospectral turbulence), the 
two integrals ( |25j ) are given by (Dung and Schlickeiser 1990, 
Schlickeiser 2002) 
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where 


G(H c ,a+,a-) =4-[a++a_+ff c (a+— a_)] (49) 

is always positive. Rp = v/|£2| = pc/(\q\Bo) is the cosmic 
ray gyroradius, Z = q/\q\ the sign of the cosmic ray particle 
and k m [ n the smallest wavenumber of the Alfven waves with 
total magnetic field strength (8 B) 2 . It is obvious that the two 
integrals © exhibit the same momentum dependence K n 
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( p/\q\ ) 2 '• Moreover, S\/Z is positive for H c > (CT+ + 

o_)/(o- — a+) and negative for H c < (c++0-)/(a_ — o+). 
Defining the maximum cosmic ray momentum with q = Qe 
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= = 1.5 • 10 16 |e|BGuG)A, 10 o —, (50) 
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where A max = 2 jck * = lOOA-ioo pc denotes the maximum 
wavelength of the Alfven waves, we obtain for the integrals 
( |47) for p < p m 
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For larger momenta p > p m the integrals (07} are zero. 
Consequently, we find for the diffusion coefficients ([33]» and 
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4. PARTICLE ACCELERATION AT RELATIVISTIC SHOCK WAVES 
We adopt the step-like shock profile ( p~6} and assume par¬ 
ticle injection S(X,p,t ) = So5(z) at the position of the shock 
only with the injection momentum spectrum S(p). Moreover, 
we assume spatially constant flow velocities and diffusion co¬ 
efficients in the upstream and downstream region. 

In the steady-state case with no losses (IRJf = 0) and a 
uniform background magnetic field (L = °°) the diffusion- 
convection transport equation (45 i in the rest frame of the 
shock wave reduces to 


d_ 

dz 


dF 

uf r — 


oco5(z) 3 


, dF p 3 F 


p 2 dp 

= S(p)8(z = 0) (54) 

We will solve Eq. ( |54} by the same method as for nonrelativis- 
tic parallel step-likeshock waves (Axford et al. 1977, Krym¬ 
sky 1987, Blandford and Ostriker 1978, Bell 1978, Drury 
1983). 

For the upstream region z > 0 the steady-state transport 
equation ([54]) yields 


Fi(z,p)=F 0 (p)exp 


Uiz 


. Hk^iJ 


(55) 


which approaches zero far upstream z —> °°. 

For the downstream region z < 0 the solution of the steady- 
state transport equation ([54]) is given by 


F 2 (z,p)=F 0 {p), 


(56) 


which is finite far downstream z —I —At the position of the 
shock 


F\(z = 0,/?) =F 2 (z = 0,p) =F 0 {p) 

the distribution function is continuous. 


(57) 


4.1. Momentum spectrum of accelerated particles at the 

shock 

The particle momentum spectrum Fo(p) at the position of 
the shock is obtained by integrating the transport equation 
( [54] ) from z = —T| to z = r| and considering the limit T| —> 0. 
THIs provides the continuity condition for the cosmic ray 
streaming density at the shock 


dF, 


rI ( U\F\ -rnfe,!-^ j -12 ^ 122 ^ 2 - 12^,2 


— T2 ( U2F2 — T2K-; 


dF 2 


ao d_ 

p 2 dp 


2 v r I dFl \ p3p{ 0 

P pzA ' ( dz J 0 + 3 


= S(p) (58) 


With the up- and downstream solutions ( [55] ) and ( [56] ) we ob¬ 
tain 


r 2 u Mp )+ u ‘ T ‘-^ ± (p + 3p2Kr! -' u ' ')FM 


*a,I 


(59) 


3 p 2 dp 
= S{p) 

where we inserted ao from Eq. (JT7T). 

According to Eqs. ( [53} we findfor the case of slab Alfven 
waves discussed in Sect. 3.4 for the ratio 


3K pz jU\ _ U\ 
Kk,i ~~ P v 


25j 


W\ 

v 


= P- 


u 1 


711 

ZR(s,a+,a-,H c ) + —- 

V 


(60) 


with the helicity dependent function 

R{s,G+,G-,H c ) — w- \ PA 
(3 — s)(5 — s) 


-C7- + H c (g+-G-)] 
(61) 

The value of the bracket of this function is not greater than 2 
for all helicity values, so the function is limited to values 


|*i, 


2(2 — s)(4 — s) 
(3 — s)(5 — s) ’ 


(62) 


which for all values of s€ [1,2) is smaller than 0.75. 

In the case of symmetric pitch-angle Fokker-Planck coeffi¬ 
cients D up {—p) = D pp (p) so that K\ = 0 the ratio ( 6 O 1 simpli¬ 
fies to 


3 K pzA U\ Uf 

—p -= 3 P^-> 

Kx.i y - 

which agrees with the limit for R — 0 of Eq. ( [60} . 
For both cases Eq. (591 then reads 


(63) 


Fo(p) 


t \fp 2 dp 
where we define the ratio 


1 d {p 2 T(p)F 0 (p)) = ® 


(64) 
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¥ = 


3t/ 2 r 2 


t/iTi -u 2 r 2 


(65) 


x J f dtt 2 S(t))ex p — y£ 


i~P dx 


and the function 


ip l J t T(x) 

so that according to Eq. © 


(76) 


+ = p(l+ Fo(p) = 


UiTi - U 2 r 2 p 2 T{p) 


'—J'dtrSWe ¥/ '^W (77) 


mcpi \A+Tp 7 me) 2 


Z/?- 


3mc(3i v/l + (p/ me) 2 


with Pi = U\/c and p = v/e. Setting 

M(p) =p 2 T(p)F Q (p) 

provides for Eq. ( |64[ > 

dM(p) \\ip 2 S(p)T(p) 

T 2 U 2 

Eq. (|68j) suggests to introduce as momentum variable 


( 66 ) 


(67) 


( 68 ) 


4.2. Monomentum injection 
For a monomomentum injection spectrum 

S{p) = S 0 8{p-po), 

the solution (|77| becomes 


Fo{p > po) = 

with the integral 


3S 0 


f/iFi -u 2 r 2 p 2 T{p) 


Po r -yt(p,po) 


(78) 


(79) 


Ap) = [ 

Jo 


dp 

o ?V) 


I(p,Po) = / 

J D 


P dx 

PO 7 » 


r 

dy 

r me 

P0 

dy 

— v 

p y 

1 + Pi \/l +y 2 [ZT? + 3pi i/l +y 2 ] 

J— V 
p y 

l + piv/l+y 2 [Zfl + 3pVl+y 2 ] 


_ fP _ t_ _ 

-Jo (1-i 2 ) [r 2 + piZ(?f + 3p2]] ’ 


implying 


dw 1 


dp T(p) 


Eq. (|68j) then becomes 


— [M(w)eV = Q( w ) 

ClXV 


with the source term 


Q(w) = 


3 p 2 ( W )S(p(w))T{p(w))eV" 
U\F\ - U 2 T 2 


Eq. (72 1 is solved for any injection spectrum S(p) by 

M(w) = / dwG{w,w)Q{w) 

Jo 

in terms of the Green s function G(w,w) obeying 
d 

dw . 

The solution of Eq. © is given by 


G{w, w )e v ' v = 8(w — w) 
© is givs 
G(w,w ) = H 


a -\|/w 


(69) 


(70) 


(71) 


(72) 


(73) 


(74) 


(75) 


= dt 


Po (l-f 2 )[f 2 + PiZ/?f + 3p 2 ]]’ 


(80) 


The integral ( [80] ) can be solved in closed form (see Appendix 
A), but asymptotic expressions for relativistic and nonrela- 
tivistic cosmic ray particle momenta can be directly deduced 
from approximating the integrand in Eq.(|8()|). 

For the differential number density of accelerated particles 
N(p) = 4 k p 2 F (p) the solution (80 1 implies 


No(p > po) = 4np 2 F(p > po) 
= 47l5 o¥ Pq v/ (p, pn ) 

u 2 r 2 t ( P ) 


(81) 


at the position of the shock z = 0, and the up- and down¬ 
stream number densities 


Ni(z > 0,p)=N 0 (p > po)exp 


Uiz 

. r l K zz,l. 


and 


N 2 {z < 0 ,p) =N 0 (p > po), 

where the ratio ( |65[ is given by 

3 


V = 


where H[x] denotes Heaviside’s step function. With Eqs. (75 i 
and (72 1 we obtain for Eq. © 


(82) 


(83) 


(84) 


i-Pf 


-1 


M(w ) =e vw [ dw Q(w ) 
Jo 


UiTi - u 2 r 2 


in terms of shock wave compression ratio r = U\/U 2 — P i /p 2 - 
In the following sections we will investigate some interest¬ 
ing limits of this solution. 
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5. SYMMETRIC PITCH ANGLE FOKKER-PLANCK COEFFICIENTS 

(R = 0 ) 

For symmetric pitch-angle Fokker-Planck coefficients (R = 
0) the integral (j80| reduces to 


I(R = 0) = [ 

■J P\ 


•p/mc 

Po/mc y 


dy 


1 +3pi(l + \) 


1 fS? ds _ 1 1 + P 2 

2 3p! + (1 +3Pj)s _ 2(! + 3Pj) n i + /| 


(85) 


with the characteristic momentum 




,, me ~ s 

l 


f V3mUi for(3i<Cl 
me for Ti 1 


( 86 ) 


Then the differential number density at the shock (81 1 be¬ 
comes the Lorentzian-type distribution function 


No(p > po) =A 0 p 




1 -P 


with 


(87) 


Aq = 


4nS 0 \\ipl 


t/ 2 r 2 (l+ 3 p2) p 2 

¥ 


P = 


2(1 + 3Pi) 


1 + — 
Pc 


+ 1 


2l p -! 

P 0 


( 88 ) 

(89) 



log(p/p c ) 

Fig. 1.— Differential number density jN} of accelerated particles at the 
shock as a function of p/p c in the case R = 0 for the adopted spectral index 
value p = 2 and injection momentum po/p c = 10 -3 . 


In Fig. 1 we illustrate the Lorentzian differential number 
density of accelerated particles at the shock as a function of 
p/p c - For particle momenta po< p< p, the Lorentzian distri¬ 
bution (87 1 increases linearly with momentum, No(po <P< 
p c ) ~ A( t p\ whereas for large momenta p > p, it approaches 
the decreasing power law distribution 


with 


No(p > p c ) ^Aop c 



(90) 


\ = 2p - 1 = 1 + 


¥ 

i+ 3 pr 


= i + 


3 

(r lv /r2^p2-l)(l + 3 p2) 


(91) 


Notice that for relativistic (Fi 1) shock speeds the charac¬ 
teristic momentum ( [86] ) coincides with me, so that in this case 
the no decreasing power law distributions for nonrelativistic 
shock accelerated particles result. 

The Lorentzian distribution ([87]) attains its maximum value 


(2p-l)P“z 

(Vo,max = Aop c ^ ' p - (92) 

at Pmax = Pc/ \/2p I. The total number of accelerated par¬ 
ticles at the shock is given by 


(Vo,tot = f JVo(p>Po) = 

■l Po 


AoPc 


l + ^) 2 


2(p — 1) L Pc . 


l-p 


4nS 0 pl _ 4 ^opg^-Pl 

u 2 r 2 cpj 

The injection rate of charged particles in a partially ion¬ 
ized (with ionization fraction x,- = n e \ /n to tt) cold upstream 
medium is related to the upstream particle number flux n t oti U\ 
as 


c $•/ \ WhotlfA . T/TeWtotlCpi s , 

Sob(p-po) = -j- °)P~P U = -2- °\P~P 0 


P\ 


P\ 


(94) 

where the injection efficiency x e <C 1 indicates the small frac¬ 
tion of ionized upstream particles being accelerated. Then the 
total number of accelerated particles at the shock (93 i is 


(Vo,tot = 47tT,T/n, ot i \jr 2 - pf 


(95) 


5.1. Nonrelativistic shock waves 

For nonrelativistic shock velocities pi <C 1, so that T| ~ 1, 
Eq. (89 1 becomes 


p~^ + l = 2r+1 

2 2(r — 1) ’ 


(96) 


where we used 


¥o = 


r — 1 


(97) 


for the limit of Eq. ( [84]) fo r pi <C 1. In this case the Lorentzian 
distribution function (|87|) reads 


No{p > po) =A 0 p 


l + ( 


V3P 


-) 2 


i me 


2r+l 
2(r—1) 


(98) 
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At momenta greater than the nonrelativistic characteristic mo¬ 
mentum p™ = \/3 p ] me this function approaches the decreas¬ 
ing power law distribution 


N,{p>p?)^A,p™(£} ** (99) 

with the spectral index 

to = (100) 

r— 1 

This spectral index agrees with the standard result for non¬ 
relativistic shocks providing > 2 for shocks in adiabatic 
electron-proton media with compression ratios r < 4 and 
qo > 3/2 for shocks in adiabatic electron-positron media with 
compression ratios r <7. 


5.2. Relativistic shock waves 

For relativistic shock velocities with pi ~ 1 and Ti 1, 
Eq. (|89|) becomes 


P- 


¥ 

8 


3 

8(r,y^3T-!)’ 


( 101 ) 


but now we have to distinguish between particle injection at 
nonrelativistic (po <C me) and at relativistic (po me) mo¬ 
menta. 

In the first case the Lorentzian distribution function m 
reads 


No{p > po) =A 0 p 


i+(— y 

me 


3 

sT^F-i-i) (102) 


At nonrelativistic particle momenta this function increases 
linearly in momentum. At relativistic particle momenta it ap¬ 
proaches the decreasing power law distribution 


6. RELATIVISTIC COSMIC RAYS 

For relativistic particle momenta the integral ( [80] ) can be 
solved for general values of the helicity dependent function 
R. We assume here that cosmic ray particles are injected at 
relativistic momenta po me with the monomomentum in¬ 
jection spectrum S(p) = So&(p — po)- With p> poi$> me the 
integral ([80]) reduces to 


I(p > po » me) ~ 


1 


l+|3iZfl + 3p 


In 


(107) 


Consequently, the solution ( |8T| ) becomes the power law dis¬ 
tribution 


No(p > po » me) ~ 


47tS 0 /W 


t/2r 2 (i + piZ/?+3Pt) \ P0 

with the power law spectral index 

¥ 




(108) 


5=1- 


= 1 + 


1 + piZ/? + 3p[ 

3 


(109) 


Fia/^-Pi-I (l +PiZ/? + 3 Pj) 


We first note that for R = 0 the power l aw so lution ( | 108] ) 
agrees with the earlie r derived expression ( |106[ ) and that the 
spectral indices (109 1 and © agree. 


6.1. Nonrelativistic shock waves 
For nonrelativistic shock velocities f/i j2 c, so that T | ~ 
1 + (Pi/2) and r 2 ~ 1 + (p//2r 2 ), the particle power law 
spectral index (109 1 to first order in pi <C 1 reduces to 


5 = 1- 


3 [1 — Z/^pr] r + 2 —3Z/?pj 


r — 1 


r — 1 


( 110 ) 


Nq (p > me) ~ Aqwc ( — ) 
V me ) 

with the spectral index 

> l) = i + 




4(TiV r 2 — 1 — 1 ’ 


(103) 


(104) 


which for Ti 1 approaches unity. We defer the discussion 
of this spectral index to the next section where the accelera¬ 
tion of relativistic cosmic rays is investigated in more detail. 

If cosmic rays are injected at relativistic momenta po 
me the power law limit ( ]90] > of the distribution function ([87 i 
holds, so that with p c ~ me again Eq. (103 i results. As po 
p c the quantity (|88]l becomes 


Aq = 


47tSo¥ 


Po 


(7 2 r 2 (l+3p 2 ) L me 

providing for the distribution (| 103]) 


2p 


(105) 


No{p > £>o » me) 


4nS Q \\tpo / p \ «*»'> 

l/ 2 r 2 (l+3p 2 ) \po) 


(106) 


To lowest order in p i we again reproduce the standard result 
for nonrelativistic shocks 


^(Pi =0) = 


r + 2 
r — 1 ’ 


(111) 


providing > 2 for shocks in adiabatic electron-proton media 
with compression ratios r <4 and q > 3/2 for shocks in adi¬ 
abatic electron-positron medi a with compression ratios r < 7. 

However, our result <| 1 10|) gives a small (for turbulence 
spectral indices s < 2) correction to this standard spectral in¬ 
dex which is different for positively (Z = 1) and negatively 
(Z = —1) charged cosmic ray particles. Depending on the 
sign of helicity dependent function R defined Eq. ( [60] ) this 
implies either a sm aller o r greater spectral index compared to 
the standard result O- With the maximum value ( [62] ) the 
correction is at most 


^ 3/? m axP1 _ 6(2 —s)(4 —s) pi 
r — 1 “ (3 — s)(5—s) r — 1’ 


( 112 ) 


which for a Kolmogorov turbulence spectral index s = 5/3 
gives 


|Afl < 3/V, xP' = 1.Q5 Pl 


r — 1 


r — 1 ’ 


(113) 
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For most adiabatic shocks with r> 1.1 this is negligibly small 
as |A^| < 10.5(3i. 


6 .2. Relativistic shock waves 


The determination of the power law spectral indices (109 1 
and © require the knowledge of the shock compression ra¬ 
tio r - P 1 /P 2 which for relativistic shocks depends for any 
given shock speed pi in a non-trivial way on the equations 
of state of the up- and downstream fluids as shown for hy- 
rodynamical shocks by Peacock (1981), Heavens and Drury 
(1988) and Kirk and Duffy (1999). The jump conditions 
for relativistic magnetohydrodynamic shocks in gyrotropic 
plasmas were studied by Double et al. (2004) and Gerbig 
and Schlickeiser (2011), including the pressure anisotropy 
X = P±/P\\ of the upstream and downstream gas pressures 


adopting adiabatic equation of states of the up- and down¬ 
stream gas with adiabatic indices K|. 2 . For a parallel ultra- 
relativistic shock (Fi > 1, pi ~ 1) Gerbig and Schlickeiser 
(2011) found the relations between the downstream parallel 
plasma beta P|i 2 = 8nPu 2 /B^, the compression ratio r = P1/P2 
and the anisotropies %i i2 . 

For illustrating our results we consider here only the case of 
an ultrarelativistic shock r 1 2S> 1 and a relativistic downstream 
medium with adiabatic index 4/3, so that pi = 3p2 (Blandford 
and McKee 1976) or r = 3. In this case we obtain for the 
power law spectral indices (| 109[> and ([9l~| 


£(n »1) = 1 
3 

i + '_ 

2 v / 2[4 + ZR]r 1 



(114) 



Fig. 2.— Power law spectral index jl 14) of relativistic particles accelerated 
at an ultrarelativistic shock for the case R = 0 as a function of the shock 
Lorentz factor I 1 . 

In Fig. 2 we calculate this spectral index for the case R = 0 
for relativistic shocks with Ti > 2, indicating spectral index 
values close to unity. Possible modifications and charge-sign 


dependencies (i.e. the case fi / 0) will be considered else¬ 
where. Due to the dominating r~ 1 dependence of q 1 the 
limit ^ — 1 is reached for Ti > 10. 

Our result of flat spectral indices with £ ~ 1 for ultra¬ 
relativistic shock disagrees strongly with the earlier estab¬ 
lished universal spectral index value q £ [2.25 — 2.30] from 
the eigenfunction and Monte Carlo simulation studies (for re¬ 
view see Kirk and Duffy 1999). As possible explanation for 
this difference we recall that our analytical solution is based 
on the two continuity conditions ( |57] > and ( [58] ) at the shock. 
These two continuity conditions are needeaas our steady- 
state diffusion-convection transport equation ( [54] ) is a second- 
order differential equation in the position coordinate z. While 
the continuity condition ( [57] ) for the particle phase density at 
the shock is also used in tKe eigenfunction solution method, 
the continuity condition ( [58] ) for the flux of particles is not 
used in that method as theFokker-Planck transport equation 
([3]) is a first-order differential equation in the position coordi¬ 
nate z. It is clear that the use of different continuity conditions 
results in different results. 

7. SUMMARY AND CONCLUSIONS 

The analytical theory of diffusive cosmic ray acceleration at 
parallel stationary shock waves with magnetostatic turbulence 
is generalized to arbitrary shock speeds V s = pic, including 
in particular relativistic speeds. This is achieved by apply¬ 
ing the diffusion approximation to the relevant Fokker-Planck 
particle transport equation formulated in the mixed comoving 
coordinate system. In this coordinate system the particle’s 
momentum coordinates p and p = p /p are taken in the rest 
frame of the streaming plasma, whereas the time and space 
coordinates are taken in the observer’s system. The Fokker- 
Planck particle transport equation contains connection coeffi¬ 
cients resulting from the coordinate transformations into this 
mixed frame which are properly included in the diffusion ap¬ 
proximation. For magnetostatic slab turbulence the diffusion- 
convection transport equation ( [45] ) for the isotropic (in the rest 
frame of the streaming plasma) part of the particle’s phase 
space density is derived for the first time for arbitrary shock 
speeds. In the limit of nonrelativistic flows the diffusion- 
convection transport equation differs from the transport equa¬ 
tion used in earlier nonrelativistic diffusive shock acceleration 
theory (Axford et al. 1977, Krymsky 1987, Blandford and 
Ostriker 1978, Bell 1978, Drury 1983) by an additional term. 
This results from our correct handling of the connection co¬ 
efficients. The additional term implies a velocity dependence 
of the acceleration rate and thus provides a major modifica¬ 
tion of the resulting differential momentum spectmm of ac¬ 
celerated particles in the nonrelativistic flow limit at nonrela¬ 
tivistic particles momenta: instead of a power law distribution 
of accelerated particles at the shock a Lorentzian distribution 
function results, which at large momenta then approaches the 
power law distribution inferred in earlier acceleration theories 
for nonrelativistic shock speeds. 

For a step-wise shock velocity profile the steady-state 
diffusion-convection transport equation is solved analytically 
again for the first time for arbitrary shock speeds follow¬ 
ing closely the solution method developed for nonrelativis¬ 
tic speeds, making use of the continuity conditions ( [57] ) and 

S for the cosmic ray phase space density and streaming 
iity at the shock. For a symmetric pitch-angle scatter¬ 
ing Fokker-Planck coefficient D^(—p) = //„,(/./) the steady- 
state solution is independent of the microphysical scattering 
details. For nonrelativistic mono-momentum particle injec- 
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tion at the shock the differential number density of acceler¬ 
ated particles is a Lorentzian-type distribution function which 
at large momenta approaches a power law distribution func¬ 
tion N(p > p c ) oc p-t= with the spectral index ^(pi) = 1 + 

[3/(Fi \Jr 2 — pj — 1)(1 +3Pj)]. For nonrelativistic (Pi <C 1) 

shock speeds this spectral index agrees with the known re¬ 
sult ^(p! <C 1) — (r + 2)/(r— 1), whereas for ultrarelativistic 
(Fj 1) shock speeds the spectral index value is close to 
unity. If particle injection occurs already at relativistic mo¬ 
menta, the steady-state solution is of power law type at all 
higher particle momenta. 

For asymmetric pitch-angle scattering Fokker-Planck co¬ 
efficient D^(—p) ^ D m (p), resulting from magnetostatic 
isospectral slab Alfven waves with non-zero values of the 
magnetic and cross helicities, the momentum spectrum of ac¬ 
celerated particles depends on the microphysical details of 
particle’s pitch angle scattering. In particular, a dependence 
of the momentum spectrum on the charge sign of the cosmic 
ray particles is found. For turbulence spectral indices s < 2, 
however, this difference is negligibly small. But steeper tur¬ 
bulence power spectra with s > 2 may provide stronger dif¬ 
ferences, which needs to be investigated in the future. Fur¬ 
ther future research topics will be concerned with the deriva¬ 
tion of time-dependent solutions of the generalized diffusion- 
convection transport equation, the study of non-stationary 
flows, the inclusion of finite up- and downstream free escape 
boundary conditions, the investigation of the influence of fi¬ 
nite frequency effects for low-frequency Alfven waves dis¬ 
carding the magnetostatic approximation, and the investiga¬ 


tion of the influence of cosmic ray momentum losses. These 
additional effects have been studied before in great detail for 
nonrelativistic shocks, and it is of high interest to investigate 
their importance for relativistic shocks. The work presented 
here provides the analytical basis for these future studies. 


I thank Msc Steffen Krakau and MSc Thorsten An- 
tecki for help with the figures and a critical reading of 
the manuscript. This work was partially supported by 
the Deutsche Forschungsgemeinschaft through grants Schl 
201/23-1 and Schl 201/29-1. 

8. APPENDIX A: THE INTEGRAL I(P. P 0 ) 

The integral ( [80] ) is be solved in closed form by 

Wi(x) 


/(P,Po)=w(|3)-w(Po), w ( x ) = 

W\(x) =2(1 — 3p[)arctan 


W 2 
ZR$\ +2x 


p,\/12 -R 2 
-f-\/12 — R 2 |"(Z/?pi - 1 - 3pj) ln(x— 1) 
— (l + 3Pi +Z/?pi) ln(l +x) 

+ (1 + 3pj) ln(3pj + Z/?p i .r E .x~ j 


W 2 = 2 


(l+3p?) 2 -fl 2 pj) \/\2-R 2 (115) 
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